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Introduction 


Efficient solution of the equations from spectral discretizations is 
essential if the high-order accuracy of these methods is to be realized in 
practice. Direct solution of these equations is rarely feasible, thus itera- 
tive techniques are required. As a consequence of the ability of spectral 
methods to capture accurately a wide bandwidth of information, the eigenvalue 
spread of a spectral operator is large, typically on the order of the square 
of a corresponding low-order finite-difference operator. Therefore, an 
explicit method would converge quite slowly, and the convergence rate would 
deteriorate very rapidly with mesh refinement. Implicit methods are not effi- 
cient, since spectral operators are full, rather than banded as in the case of 
finite difference. 

An alternative to pure explicit iterative schemes is a preconditioned 
scheme, where explicit iteration of some form is driven not by the spectral 
residual, but by the residual obtained after some processing is applied to 
reduce its eigenvalue spread. As will be shown later, this process can also 
be thought of as an approximate-implicit scheme, which can give some Insight 
into a relevant preconditioning operator. 
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Preconditioned-iteration schemes for spectral collocation discretization 
of second-order equations are well-known and proven. For instance, the 
time-accurate incompressible Navier-Stokes simulations of certain fluid 
mechanical phenomena, in which Chebyshev collocation is used in two coordinate 
directions, require the solution of a number of Helmholtz or Poisson equations 
per time step (ref. 1). The preconditioned scheme used in that work requires 
less than one second to accomplish such a solution on a 65 x 65 mesh on a 
CYBER 205 vector computer. The preconditioning operator for these 
second-order equations is the low-order central finite difference operator, 
using the Chebyshev collocation points as its mesh. Orszag (ref. 2), 
originally suggested such preconditioning and provided some analysis for the 
case of Fourier discretization. 

The steady -state compressible Euler or Navier-Stokes equations at high 
Reynolds numbers, however, are advecti on -dominated. A spectral solution 
technique for such equations must therefore deal with operators which are 
predominantly first-order. Since inversion of the finite -difference 
preconditioning operator is (as will be shown) related to the 
finite-difference solution to the original problem, one would like to draw on 
finite-difference experience in solving compressible Navier-Stokes problems to 
yield efficient solution of the preconditioning operator. Variations of the 
Beam and Warming scheme (refs. 3, 4) are popular for such solutions. However, 
the advecti on terms in this scheme are central differenced. Elementary 
analysis of preconditioning first-order Fourier discretization with central 
finite -difference on the collocation mesh indicate that the convergence rate 
of such a scheme would be unacceptable; the eigenvalue spread of this 
preconditioned operator is unbounded. A similar situation appears in the case 
of Chebyshev discretization, as will be shown here. 
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A preconditioning scheme for first-order Chebyshev collocation operators 
is proposed herein, in which the central finite-difference mesh is finer than 
the collocation mesh. (The authors later found that Orszag conjectured such a 
preconditioner in an early paper for Fourier series, although no rigorous 
analysis was given (ref. 2).) Details of the proper techniques for 
transferring information between the meshes is given here, and the scheme is 
analyzed by examination of the eigenvalue spectra of the preconditioned 
operators, corresponding to a pure first-order and to an advecti on-dominated 
advection/dif fusion problem, both with realistic boundary conditions. The 
effect of artificial viscosity required in the inversion of the 
finite-difference operator is examined. A second preconditioning scheme, 
involving a high -order upwind finite -difference operator of the van Leer type 
(ref. 5), is also analyzed to provide a comparison with the present scheme. 
Finally, the performance of the present scheme is verified by application to 
several test problems. 


Overview of Preconditioning Iteration Schemes 


Consider the following linear equation 

L sp u ■ f (1) 

where the operator L sp is derived from spectral collocation of a 
differential equation. An iterative scheme is to be used to solve this 
equation. Given a current estimate of the solution u n at iterate "n", a 
simple Richardson iteration scheme for computing a better estimate u n+1 is 

u" +1 = u n - u (L u n - f ) (2) 

sp 

where the scalar relaxation factor w may be chosen either experimentally or 
via a requirement that some norm of the residual 
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be minimized at each step. Rewriting the scheme (2) as 

u n+1 - u n = Au n = -W R n (4) 

shows the explicit nature of this iteration. Preconditioning involves 
choosing an operator M which is more easily invertible than L sp , and for 
which the scheme 

A u n = -w M' 1 R n (5) 

converges more rapidly than the scheme (4). The convergence rate of such 

schemes is quantified in the following way: expand Eq. (5) as 


u n+1 = (I - a) M" 1 L ) u n - as M" 1 f (6) 

sp 

where I is the identity operator. Subtracting the discrete solution to Eq. 
(1) (the desired solution u) from both sides of (6), and adding to the right 
hand side 


o> M" 1 (L c u - f) , (7) 

bp 

which is equal to zero by (1), yields an equation for the discrete error: 

(u n+1 - u) = (I - to M" 1 L gp ) (u n - u) (8) 

For the preconditioner scheme (5) to be convergent, all norms of (u n - u) 
must decrease; thus 


ii I - oi M" 1 L I < 1 (9) 

sp 

is required. Given a preconditioning operator M, the relaxation factor w is 
used to satisfy Eq. (9). 


* 
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From Eq. (9) may be seen what constitutes an effective preconditioner: 
if the eigenvalues of M" 1 L S p are clustered in a unit circle centered at 1 
in the complex plane, then all error components will converge at nearly the 
same rate, and an optimal w may be chosen which will yield rapid convergence. 

A slow scheme is characterized by a wide spread of eigenvalues from such a 
clustered pattern. 

As an aside, a preconditioned iteration scheme may be looked at as an 
approximate implicit scheme in the following way. Ideally, given an estimate 
u", one wants the residual at the next iterate to be zero; expanding 

Rn+1 = L $p (u n + Au 11 ) - f = 0 (10) 

or 

L sp iu = -R" (11) 

However, L sp is difficult to invert; approximate it on the left by a more 
easily-inverted operator M, giving 

Au n = -0) M" 1 R n (12) 

which is identical to Eq. (5). Thus the better aT* M approximates L sp , 
the faster the scheme will converge. 

The inverse of the preconditioning operator may be obtained and used 
directly, as implied by Eq. (5) or the preconditioning equation 

M au 11 = -id R n (13) 

may be inverted iteratively. This procedure is beneficial when the spectral 
operator L sp is nonlinear, necessitating a different operator M at each 
iterate, or when M is still expensive to invert. An example of the latter 
case is when Eq. (1) results from spectral discretization of a PDE in two or 
more dimensions. The corresponding finite-difference discretization operator 
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works well as a preconditioner (with the modifications described in the next 
section for first-order operators), but may still be expensive to invert, due 
to its size. An iterative technique for computing the Au n is more efficient 
than direct inversion for this case, using, for instance, an approximate-fac- 
torization scheme. Experience indicates that the preconditioner-inversion 
iteration need not be driven to convergence for the overall preconditioned 
scheme (Eq. (5)) to converge rapidly. 

Preconditioners for First-Order Spectral Operators 

As mentioned in the introduction, the preconditioning scheme investigated 
here uses a central finite-difference operator, constructed on a mesh finer 
than that of the spectral discretization. A simple example will show why pre- 
conditioning using central finite-differences on the collocation mesh is inad- 
equate. For the model scalar problem U x = f with periodic boundary condi- 
tions on [0, 2tt ) , the eigenvalues of the Fourier collocation operator 3/3x are 
iaAx, where a is the wavenumber and Ax is the constant collocation mesh 
spacing. The product aAx falls in the range 0 <|aAx|< ir. The corresponding 
eigenvalues for the central -difference operator are i»sin (aAx). Note that as 
aAx-Mr the ratio of these eigenvalues is unbounded. This ratio corresponds to 
the eigenvalues of the preconditioned operator, denoted M"* L S p in the 
previous section. Such a preconditioned scheme is thus unconditionally 
unstable, with unbounded growth of the highest wavenumber error components. 

The use of a finite preconditioning grid averts this unbounded component by 
introducing some natural dissipation in the preconditioner at the highest wave 
number of the spectral operator. 

Because of the difference between the spectral and preconditioner meshes. 
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a scheme is needed for transferring information between them. In multi grid 
terminology, a prolongation operator is needed to transfer the spectral 
residual which acts as the source term of the preconditioning equation (Eq. 
(13)), and a restriction operator is needed to compute the solution updates on 
the spectral mesh. Naturally, it is desirable to transfer as much information 
as possible form the spectral residual to the preconditioning equation; 
spectral interpolation is therefore used for the prolongation operator. In 
the restriction operation, however, aliasing of correction components, with 
wavenumber higher than that of the spectral mesh must be avoided. Therefore, 
spectral restriction cannot be used; low-order Lagrange interpolation is used 
here. 

The preconditioning scheme proposed here proceeds as follows for the 
spectral discretization of a simple model problem U x = f. At each iterate, 
compute the residual (Eq. (3)) defined on N sp points in the domain. This 
information is transferred to the (finer) preconditioning mesh via the 

spectral interpolation operator I^j . Denoting 

R - I™ R (14) 

the preconditioning equation becomes 

M Au = R (15) 
where au is the update on the preconditioning mesh. For the model problem 
U x = f in x e [-1,1] considered here, the preconditioner M must approximate 
the first derivative operator via central finite differences. We use 

« = 8 X " - e 6 X ° X (16) 

using standard divided-difference notation. The second-difference term is 
required to avert the odd-even uncoupling of the pure central first-difference 
operator, and to aid in its inversion. Some type of artificial viscosity is 
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essentially always required when inverting this type of operator, including 
application in the Beam and Warming scheme. Note that this artificial 
viscosity does not effect the spectral solution, being confined to the precon- 
ditioner. The magnitude of the parameter is quite small, generally 10~ 3 or 
less; its effect will be discussed in the next section. 

Since we are concerned here with general, nonperiodic problems, the 
effect of boundary conditions on the preconditioning scheme is relevant. The 
appropriate spectral boundary condition for this problem is Dirichlet on one 
end of the domain (x = -1). The other end (x = 1) requires no boundary condi- 
tion; the spectral scheme is used to enforce the equation at that point. On 
the other hand, extrapolation conditions are commonly used in finite-differ- 
ence discretizations of such problems. Such a boundary condition is used in 
the preconditioner as analyzed in the next section. 

The final step of this preconditioning scheme is to carry the update 
information to the spectral mesh. The iterate thus concludes with 

au = Ip^ aU 
and 

u n+1 = U 11 + 0 ) au (17) 

where the operator 1^ uses low-order Lagrange interpolation. For analysis 

purposes, the above sequence may be collected into a single operator: 

H- 1 = 7*j> H- 1 I™ (18) 

This is, in effect, the preconditioning operator applied in the above scheme. 

In practice on a real problem of interest, the finite-difference operator 
M is also too large to invert directly at reasonable cost. Preconditioning 
operators are therefore inverted iteratively in all but the simplest one- 
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dimensional test problems. Experience with both first- and second-order 
operators indicates that full machine-zero convergence of this iterative 
inversion is not necessary; depending on the nature of the iterative scheme, 
only two orders of magnitude reduction in the residual is necessary. 

One preconditioning operator which is attractive from the standpoint of 
there being a rapid and efficient inversion technique involves high-order 
upwind differencing of the van Leer type (ref. 5). In this scheme, the flux 
is split according to direction of propagation, and appropriate upwind 
differences are taken. The inherent artificial viscosity in such schemes has 
been shown to provide enhanced robustness, and allows for nearly 
fully-implicit treatment of the artificial viscosity during iteration for even 
complex problems, which speeds convergence. 

A simple first derivative, taken as a forward-propagating flux, is 
approximated as 

6 u — ~ u i-i/2^ 1 ^ (19) 

where 

u| ± l /2 = u. ± */4 [(1 + k) V + (1 ± k) A] u (20) 

The operators v and a are undivided backward and forward differences, 
respectively. The parameter k determines the accuracy of the operator: k=~l 

corresponds to a fully-upwind scheme, k=l to a central difference operator, 
k=l/3 to a third-order upwind biased scheme. 
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Results 


To investigate the effectiveness of the proposed preconditioning scheme, 
the eigenvalues of the operator M _1 L sp were computed for a number of 
combinations of N sp and Np^, for the model problem 

U = sin trx, x e [-1,1] (21) 

A 

U(-l) = 0 

These results were compared qualitatively with convergence rate observations 
from actual implementation of the scheme. 

Some basic characteristics of the eigenvalue spectra for this problem are 
as follows: The eigenvalues of the spectral operator L sp alone are 

O 

dominantly imaginary, with the magnitude of the largest growing as N^ p 

When the spectral operator is preconditioned with the central finite- 
difference operator on the same mesh, the eigenvalues of the overall operator 
become dominantly real, which is a desirable feature for use in an iterative 
scheme. The real parts of these eigenvalues are all positive, the smallest 
being near 1 and the largest on the order of 50 for the grids investigated. 
The eigenvector associated with this largest eigenvalue is highly oscillatory, 
as expected from the discussion in the previous section. As Npp is 
increased relative to N sp , this eigenvalue pattern generally collapses onto 
the point (1,0) in the complex plane, with both real and imaginary parts 
decreasing. This collapse is at first rapid as Npp is increased from N sp , 
then slows as the eigenvalue with the largest real part moves to the interior 
of the unit circle centered at 1. There is a single exception to this 
clustered pattern, that of a small (0( . 1 ) ) eigenvalue which remain essentially 
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fixed for all Npo < N sp . Its eigenvector indicates that this eigenvalue 
is associated with the different conditions used at x = 1 by M and L sp . In 
general, though, the eigenvalues of the preconditioned operator are strongly 
clustered, thus one would expect rapid convergence of the iterative scheme. 

In Table I are shown the maximum and minimum eigenvalues for the precon- 
ditioned operator from the test problem described above, for various combina- 
tions of N sp and Npo. The range of N sp considered covers what we expect 
to be required in spectral discretizations of practical aerodynamic problems. 
As can be seen, the maximum eigenvalue is large when Npp = N sp , and drops 
rapidly as Npo is increased. The operator appears well -conditioned, that 
is, the eigenvalues are tightly grouped inside the unit circle centered at 1, 
when Npg 1.5 N sp for all of the grids considered. These results were 
produced with first-order Lagrange interpolation in the restriction operator, 
and with the artificial -viscosity coefficient used in M fixed at 10“3. 
Virtually identical results were obtained for e = 10-4, anc i f or second-order 
restriction. 

In Table II are shown the convergence rates observed from application of 
this preconditioning scheme to the model problem Eq. (21). Tabulated are the 
average reduction per iteration of both maximum residual and maximum error, 
taken over a reduction in residual of at least eleven orders of magnitude. 
Results from two iterative schemes are shown in Table II: the first using a 
fixed (o, the near-optimum of which was chosen by trial -and-error; the second 
in which the w was chosen by the requirement that the Lz-norm of the resi- 
dual be minimized at each iteration, given an update "direction." This latter 
minimum-residual scheme has been found to generally coincide well with the 
results of Table I, with fixed-t) values near 1 when all eigenvalues lie inside 
the aformentioned circle; reduced values of u> are required or non -convergence 
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is seen when conditioning of the operator is inadequate. This test also 
indicates that Npg > 1.5 N sp yields good results, with convergence rate 
increasing as Npg is increased beyond this point. Note also that the con- 
vergence rates for the N sp = 45, NpQ =75 case, where Npp = 5/3 N sp , 
is faster than the case N sp = 33, Npp = 55, where the ratio is again 5/3. 

<• 

This indicates that the convergence rate of the scheme does not decrease with 
spectral mesh refinement, which is a very desirable property. 

A similar eigenvalue spectrum behavior is seen in Table III for an appli- 
cation of the present preconditioner to the advecti on-diffusion problem 

U x - v U xx = f , x e [-1,1] (22) 

U(-l) = a , U(l) = b 

The preconditioning operator used here is 

M = 6° - (e + v) 6° x (23) 

The results presented in the first part of Table III are for e = v = 10~ 3 ; 
that is, some additional, artificial viscosity is applied in the 
preconditioning operator augmenting the "natural" viscosity. Again, the 
preconditioner effectively clusters the eigenvalues of the operator. In the 
latter part of Table III are shown the eigenvalues for a case with the same 
physical viscosity, v = 10~ 3 , but with double the artificial viscosity of 
the former case; e = 2 x lO"' 3 . Only a small effect is seen once 
Npg > 1.5 N sp . 

Maximum and minimum eigenvalues for first-, second- and third-order up- 
wind schemes preconditioning Eq. (22) with v = 10~ 3 are shown in Table IV. * 

Note that although the eigenvalues are well-behaved in magnitude, the imagin- 
ary parts are considerably larger than in the previous scheme. Convergence- 
acceleration methods such as mi nimumal -residual choice of w are less effective 
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on operators with eigenvalues having large imaginary parts, or may be unstable 
or oscillatory during convergence. Such a feature is therefore not good in a 
precondi ti oni ng appl i cati on . 

To illustrate the differences in eigenvalue spectra between the present, 
central ^differences preconditioning and the high*»order upwind preconditioning, 
eigenvalues are plotted in Figs. 1^6 for Eq. 22, v = 1(T 5 . The artificial 
viscosity for the central ^difference scheme was e = 1(T 4 . Shown for 
N S p = 45 and 60 are the spectra (positive imaginary parts) without 
preconditioning (Figs. 1 and 4), for NpQ = 1.2 N sp (Fig. 2a, 5a), for 
Npo = 1.5 N S p (Figs. 2b, 5b), and for first**, second**, and third**order 
upwind (Figs. 3, 6). Note that except for a few spurious eigenvalues, the 
spectra for central **difference preconditioning are well ^contained and have 
small imaginary parts, whereas those of the upwind preconditioners have large 
imaginary parts, increasing as the real parts tend to zero. As stated before, 
such behavior is detrimental to iterative convergence of the overall 
preconditioned scheme. 


Conclusions 

We have demonstrated an effective preconditioner for the operators which 
result from collocation discretization of first derivatives. Results of an 
eigenvalue analysis of the precondi tionedMterati on scheme agreed well with 
actual results for a simple test problem. The preconditioner performed well, 
producing rapid convergence. 

The proposed scheme is advantageous in application as a preconditioner 
for spectral collocation solutions of the compressible Navier^Stokes equations 
for several reasons. First, since solution of the preconditioning equation is 
related to solution of the original equation via central finite differences. 
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one can draw on well -developed techniques, such as the Beam and Warming 
scheme. Second, physical ly-relevant boundary conditions can be imposed in the 
preconditioner, so that boundary points are included consistently in the 
update scheme. For outflow boundaries, when the spectral scheme merely 
enforces the given equation, it was found that simple extrapolation conditions 
in the preconditioner were adequate. Finally, it was found that the conver- 
gence rate of the preconditioned scheme was essentially independent of the 
spectral mesh size, when the preconditioning mesh was fine enough. The degree 
of refinement over the spectral mesh required in the preconditioner was not 
large; 50% more points was found to be adequate. Such a mesh would still be 
considered quite coarse from finite -difference standards, so the greatly 
increased accuracy of the spectral method can be economically realized. It 
was also shown that high-order upwind schemes, despite a potential advantage 
in more efficient inversion of the preconditioning operation, probably will 
not perform as well, due to the large imaginary parts of the eigenvalues of 
the preconditioned operator. 

Currently underway are investigations using two-dimensional test problems 
to verify in that setting the findings presented here. One aerodynamic appli- 
cation planned is that of supersonic flow over a blunt leading edge, using 
bow-shock fitting in a compressible Navier-Stokes formulation. This problem 
has been investigated using a finite-difference Beam and Warming code (ref. 
4); it was found that truncation error compromised heat-transfer predictions, 
thus a higher-order accuracy technique is required. Other applications 
planned involve basic studies of fluid-dynamic stability of compressible flow, 
such as high-speed boundary layers. Spectral methods are needed for such sim- 
ulations due to the stringent accuracy requirements of wave propagation pre- 
diction in this setting. 
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Table I. Maximum and minimum eigenvalues for various N sp , Np D 
preconditioning Eq. (21). (Real, imaginary parts) 


N sp 

N F d 

x max 

x min 

9 

9 

21.48, 0.0 

1., 0.0 


14 

1.208, ± .046 

.339, ± .310 


18 

1.054, ± .018 

.173, 0.0 

17 

17 

47.46 , 0.0 

1., 0.0 


21 

2.333,0.0 

.276, 0.0 


26 

1.364,0.0 

.218, 0.0 


34 

1.136,0.0 

.082, 0.0 

33 

33 

56.28, 0.0 

1., 0.0 


40 

2.178, 0.0 

.097, 0.0 


55 

1.365, ± .013 

.056, 0.0 


66 

1.348, 0.0 

.049, 0.0 

45 

45 

37.42, 0.0 

1., ± .002 


54 

2.456, 0.0 

.073, 0.0 


68 

1.53, 0.0 

.049, 0.0 


75 

1.390, ± .062 

.045, 0.0 
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Table II. Convergence rates of preconditioned schemes for Eq. (21). 
p(R), p (a ) are average reduction per iteration of residual, error, 
respectively. 

fixed (o MR 


^sp 

n FD 

U) 

p(R) 

p(A) 

p(R) 

P ( A ) 

9 

9 


=41^ 


Hi 5 * 

Hi 5 * 


14 

.4 

.616 

.601 

.508 

.501 


18 

.5 

.496 

.501 

.391 

.383 

17 

17 




^N 5 * 

HI 5 * 


21 

.5 

.828 

.834 

HI 5 * 

HI 5 * 


26 

.4 

.639 

.615 

.225 

.257 


34 

.5 

.515 

.508 



33 

33 


HI 5 * 


HI 5 * 

HI 5 * 


40 

.5 

.609 

.576 

HI 5 * 

HI 5 * 


55 

1.0 

.419 

.386 

HI 5 * 

HI 5 * 


66 

.9 

.267 

.207 

.120 

.112 

45 

45 


Hi 5 * 


s-N 5 * 

^N 5 * 


54 

.5 

.619 

.566 

HI 5 * 

HI 5 * 


68 

1.0 

.553 

.494 

.283 

.224 


75 

1.0 

.337 

.291 

.239 

.196 


* 


j 
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Table III. Maximum and minimum eigenvalues for various N sp , Npp 
preconditioning Eq. (22). (Real, Imaginary parts) 



N sp 

n FD 

x max 


x min 


33 

33 

29.152, 

0.0 

.753, 

0.0 


40 

2.190, ± 

1.269 

.103, 

0.0 


55 

1.360, ± 

.056 

.067, 

0.0 


66 

1.169, ± 

.120 

.062, 

0.0 

45 

45 

13.024, ± 

7.955 

.754, 

0.0 


54 

2.329, ± 

.140 

.124, 

0.0 


68 

1.502, ± 

.092 

.086, 

0.0 


75 

1.381, ± 

.079 

.077, 

0.0 


v = lO" 3 , 

e = 2 x 

10" 3 




Nsp 

n fd 

x max 

x min 

33 

33 

21.590, 0.0 

.997, 0.0 


40 

2.316, 0.0 

.072, 0.0 


55 

1.400, 0.0 

.070, 0.0 


66 

1.261, 0.0 

.046, 0.0 

45 

45 

109.8, 0.0 

.998, 0.0 


54 

2.709, ± 1.058 

.119, 0.0 


68 

1.591, 0.0 

.082, 0.0 


75 

1.454, 0.0 

.074, 0.0 
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Table IV. Maximum and minimum eigenvalues for upwind schemes preconditioning 
Eq. (22) (Real, imaginary parts) 


v = 1CT 3 


N 

order 

x max 

x min 

33 

1 

.231, ± 1.359 

.251, ± .448 


2 

1.000, 0.0 

.173, ± .258 


3 

.651, ± 1.843 

±0.0 

45 

1 

.211, ± 1.386 

.416, ± .435 


2 

1.000, 0.0 

.268, ± .285 


3 

.547, ± 1.927 

.999, 0.0 


a 
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° NS= 60 



Fig. 4. Eigenvalues of spectral discretization of Eq . 22, 
v = 10~ 5 . 


° NS= 60. NF= 72 



(a) N fd - 1.2 Ngp 

° NS= 60. NF= 90 



(b) N F q - 1.5 Ngp 

Fig. 5. Eigenvalues of spectral discretization of Eq. 22, 
v = 10- 5 , preconditioned with present scheme, 
e = 10 -4 . 
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° N= 60. 1st 



(a) first order 


o N= 60. 2nd 



real 


(b) second order 


° N= 60. 3rd 


O 



(c) third order 


Fig. 6. Eigenvalues of spectral discretization of Eq. 22, 
v = 10 -5 , preconditioned with upwind scheme. 
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